// Copyright (c) 2011 Bo Zhou<Bo.Schwarzstein@gmail.com>
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//		http://www.apache.org/licenses/LICENSE-2.0 
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

/**
 * Transform a vector by 3 basis vectors.
 */
vector local_spherical(varying vector v, x, y, z)
{
	float a = v[0] * x[0] + v[1] * y[0] + v[2] * z[0];
	float b = v[0] * x[1] + v[1] * y[1] + v[2] * z[1];
	float c = v[0] * x[2] + v[1] * y[2] + v[2] * z[2];

	return vector(atan(b, a), PI * 0.5 - acos(c), 0);
}
 
/**
 * Unit-height zero-mean Gaussian function.
 */
float gaussian(float x, mu, sigma)
{
	return exp(- pow(x - mu, 2.0) * 0.5 / pow(sigma, 2.0));
}

/**
 * Jitted vector.
 */
vector jitted(vector vec; float a, b, c)
{
	float amp = a * sin(2 * PI * b * v + c);
	point Po = transform("object", P);
	normal jitted = noise(Po);
	return normalize(vec + normalize(jitted) * amp);
}

color AFMarschner
(
	float phi, thetaD, thetaH;

	float intensityR;
	color colorR;
	float longitudinalShiftR;
	float longitudinalWidthR;
	
	float intensityTT;
	color colorTT;
	float longitudinalShiftTT;
	float longitudinalWidthTT;
	float azimuthalWidthTT;
	
	float intensityTRT;
	color colorTRT;
	float longitudinalShiftTRT;
	float longitudinalWidthTRT;
	
	float intensityG;
	float azimuthalShiftG;
	float azimuthalWidthG;
	
	output varying color R;
	output varying color TT;
	output varying color TRT;
)
{
	// R
	//
	float M_R = gaussian(thetaH, radians(longitudinalShiftR), radians(longitudinalWidthR));
	float N_R = cos(phi * 0.5);
	
	R = intensityR * colorR * M_R * N_R;
	
	// TT
	//
	float M_TT = gaussian(thetaH, radians(longitudinalShiftTT), radians(longitudinalWidthTT));
	float N_TT = gaussian(PI, abs(phi), radians(azimuthalWidthTT));
	
	TT = intensityTT * colorTT * M_TT * N_TT;
	
	// TRT
	//
	float M_TRT = gaussian(thetaH, radians(longitudinalShiftTRT), radians(longitudinalWidthTRT));
	float N_TRGNG = cos(phi * 0.5);
	float N_G = intensityG * gaussian(radians(azimuthalShiftG), abs(phi), radians(azimuthalWidthG));
	float N_TRT = N_TRGNG + N_G;
	
	TRT = intensityTRT * colorTRT * M_TRT * N_TRT;
	
	// Normalize.
	//
	float cosThetaD = cos(thetaD);
	float scale = 1.0 / pow(cosThetaD, 2);

	return (R + TT + TRT) * scale;
}